Практическое задание №4

Плискин Александр

In [1]:
import numpy
import time
import scipy.special
import scipy.misc
import npy2cube

Параметры

n - основное число (1,2,3..)
l - орбитальное число (0,1,2.. n-1)
m - магнитное число (-l..+l)

Переход от декартовых координат в сферические

$ r (x,y,z)=\sqrt{x^2+y^2+z^2} $ - расстояние до начала координат.

$ \theta(x,y,z)=arccos\Big(\frac{z}{r(x,y,z)}\Big) $ - зенитный угол.

$ \phi(x,y,z)=arctg\Big(\frac{y}{x}\Big) $ - азимутальный угол.

Волновая функция:

$ a_0 $ - уменьшенный Боровский радиус.

$ R(r,n,\ell) = = \Big(\frac{2r}{na_0}\Big)^\ell e^{-\frac{r}{na_0}}L^{2\ell+1}_{n-\ell-1}\frac{2r}{na_0} $ - радиальная часть функции, где $ L_{n-\ell-1}^{2\ell+1}(...) $ - обобщённый полином Лагерра степени $ n-\ell-1 $.

$ WF(r,\theta,\phi,n,\ell,m) = R\cdot Y^m_\ell(\theta,\phi) $ - Волновая функция, заданая как произведение радиальной и угловой частей, где $ Y_{\ell}^{m}(\theta, \phi) $ - сферическая гармоника степени $ \ell $ порядка $ m $.

$ absWF(r,\theta,\phi,n,\ell,m) = |WF(n,\ell,m,r,\theta,\phi)|^2 $ - плотность вероятности нахождения частицы в данной точке пространства в данный момент времени.

In [2]:
def w(n,l,m,d):

    x,y,z = numpy.mgrid[-d:d:30j,-d:d:30j,-d:d:30j]

    r = lambda x,y,z: numpy.sqrt(x**2+y**2+z**2)
    theta = lambda x,y,z: numpy.arccos(z/r(x,y,z))
    phi = lambda x,y,z: numpy.arctan(y/x)

    a0 = 1.

    R = lambda r,n,l: (2*r/n/a0)**l * numpy.exp(-r/n/a0) * scipy.special.genlaguerre(n-l-1,2*l+1)(2*r/n/a0)
    WF = lambda r,theta,phi,n,l,m: R(r,n,l) * scipy.special.sph_harm(m,l,phi,theta)
    absWF = lambda r,theta,phi,n,l,m: numpy.absolute(WF(r,theta,phi,n,l,m))**2

    return WF(r(x,y,z),theta(x,y,z),phi(x,y,z),n,l,m)
Расчет значений для первых трех уровней. Создание Cube-файлов.

Также пришлось немного изменить скрипт npy2cude для совместимости python3

In [3]:
d = 30
step = float(2.*d/29)

# Зададим цикл по перебору квантовых чисел
for n in range(0,4):
    for l in range(0,n):
        for m in range(0,l+1,1):
            grid = w(n,l,m,d)
            name ='%s-%s-%s' % (n,l,m)
            npy2cube.npy2cube(grid,(-d,-d,-d),(step,step,step),name+'.cube')
/Users/alex/Desktop/ws/hse/golovin/task4/npy2cube.py:41: ComplexWarning: Casting complex values to real discards the imaginary part
  cube_file.write("%f " %(float(grid[x, y, z])))
/Users/alex/Desktop/ws/hse/golovin/task4/npy2cube.py:44: ComplexWarning: Casting complex values to real discards the imaginary part
  cube_file.write("%f\n" %(float(grid[x, y, z])))

PyMOL

In [4]:
from IPython.display import display,Image
from xmlrpc.client import ServerProxy
cmd = ServerProxy(uri='http://localhost:9124/RPC2')
In [5]:
# Откуда эти цифры? 
# vol obj -> coloring -> panel -> *manual change* -> get colors as script
# cmd.volume_ramp_new('name of profile', [\
#     pointi:  thickness,   R,   G,   B,   opacity,\ ...]) 

cmd.reinitialize()
cmd.volume_ramp_new('ramp007', [\
        0.002, 0.00, 0.00, 1.00, 0.00, \
        0.01,  0.00, 1.00, 1.00, 0.20, \
        0.015, 0.00, 0.00, 1.00, 0.00, \
    ])
cmd.set('ray_volume')
cmd.load('/Users/alex/Desktop/ws/hse/golovin/task4/1-0-0.cube')
cmd.volume('vol', '1-0-0')
cmd.volume_color ('vol','ramp007')
cmd.ray(300, 300)
cmd.png('/Users/alex/Desktop/ws/hse/golovin/task4/1-0-0.png')
time.sleep(1)
Image('/Users/alex/Desktop/ws/hse/golovin/task4/1-0-0.png')
Out[5]:
In [6]:
cmd.reinitialize()
time.sleep(2)
cmd.set('ray_volume')
cmd.volume_ramp_new('ramp007', [\
        -0.05, 0.00, 0.00, 1.00, 0.00, \
        -0.01,  0.00, 1.00, 1.00, 0.20, \
        -0.015, 0.00, 0.00, 1.00, 0.00, \
        0.005, 0.00, 0.00, 1.00, 0.00, \
        0.01,  0.00, 1.00, 1.00, 0.20, \
        0.01, 0.00, 0.00, 1.00, 0.00, \
        ])
names = []
for n in range(0,4):
    for l in range(0,n):
        for m in range(0,l+1,1):
            name='%s-%s-%s' % (n,l,m)
            file = '/Users/alex/Desktop/ws/hse/golovin/task4/' + name + '.cube'
            cmd.load(file)
            cmd.volume('vol', name)
            cmd.volume_color ('vol','ramp007')
            time.sleep(10)
            cmd.png('/Users/alex/Desktop/ws/hse/golovin/task4/' + name + '.png')
            time.sleep(5)
            names.append(name)
In [7]:
for name in names:
    display(Image('/Users/alex/Desktop/ws/hse/golovin/task4/' + name + '.png'))

ORCA

определенная схожесть имеется

In [8]:
cmd.set('ray_volume')
cmd.volume_ramp_new('ramp007', [\
        -0.05, 0.00, 0.00, 1.00, 0.00, \
        -0.01,  0.00, 1.00, 1.00, 0.20, \
        -0.015, 0.00, 0.00, 1.00, 0.00, \
        0.005, 0.00, 0.00, 1.00, 0.00, \
        0.01,  0.00, 1.00, 1.00, 0.20, \
        0.015, 0.00, 0.00, 1.00, 0.00, \
        ])
for i in range(1,5):
    name='H-%s'% i
    file = '/Users/alex/Desktop/ws/hse/golovin/task4/' + name + '.cube'  
    cmd.load(file)
    cmd.volume('vol', name)
    time.sleep(10)
    cmd.png('/Users/alex/Desktop/ws/hse/golovin/task4/' + name + '.png')
    time.sleep(5)
In [9]:
for n in range(1,5):
    name = 'H-%s' % n
    display(Image('/Users/alex/Desktop/ws/hse/golovin/task4/' + name + '.png'), name)
'H-1'
'H-2'
'H-3'
'H-4'